//============================================================================
// Name        : csiena
// Author      : Josh Lospinoso
// Copyright   : MIT License
// Description : Austere, limited implementation of the R package RSiena in
//				 C++. Used for fast maximum likelihood and method of moments
//				 estimation of stochastic actor oriented models (SAOMs)
//				 without any frills or user experience niceties.
//============================================================================

#include "chain.h"

namespace csiena {

Chain::Chain(int actors) {
	this->linkEvaluationFunction = new LinkEvaluation(actors);
	this->crimpEvaluationFunction = new CrimpEvaluation(actors);
	this->startingDigraph = new Digraph(actors);
	this->currentDigraph = new Digraph(actors);
	this->currentElement = 0;

	this->betaScore = new double[BETA_LENGTH];
	this->betaInformation = new double[BETA_LENGTH * BETA_LENGTH];
	#if GAMMA_LENGTH > 0
	this->gammaScore = new double[GAMMA_LENGTH];
	this->gammaInformation = new double[GAMMA_LENGTH*GAMMA_LENGTH];
	#endif
	this->actors = actors;
	this->currentElementIndex = 0;
	this->chainSize = 0;
	this->diagonalElements=0;
	this->crimps = 0;
	this->alpha = 1;
}

Chain::Chain(const Chain &donorChain) {
	// Basic memory allocation
	this->betaScore= new double[BETA_LENGTH];
	this->betaInformation=new double[BETA_LENGTH * BETA_LENGTH];
	#if GAMMA_LENGTH > 0
	this->gammaScore = new double[GAMMA_LENGTH];
	this->gammaInformation = new double[GAMMA_LENGTH*GAMMA_LENGTH];
	#endif
	this->actors = donorChain.actors;
	this->logLikelihood = donorChain.logLikelihood;
	for (int i=0; i<BETA_LENGTH; i++) {
		this->betaScore[i] = donorChain.betaScore[i];
		for (int j=0; j<BETA_LENGTH; j++) {
			this->betaInformation[i*BETA_LENGTH+j] = donorChain.betaInformation[i*BETA_LENGTH+j];
		}
	}
	for (int i=0; i<GAMMA_LENGTH; i++) {
		for (int j=0; j<GAMMA_LENGTH; j++) {
			this->gammaInformation[i*GAMMA_LENGTH+j] = donorChain.gammaInformation[i*GAMMA_LENGTH+j];
		}
	}
	this->chainSize = donorChain.chainSize;
	this->diagonalElements = donorChain.diagonalElements;
	this->alpha = donorChain.alpha;

	//Need to be copied via copy constructors
	this->linkEvaluationFunction = new LinkEvaluation(*(donorChain.linkEvaluationFunction));
	this->crimpEvaluationFunction = new CrimpEvaluation(*(donorChain.crimpEvaluationFunction));
	this->startingDigraph = new Digraph(*(donorChain.startingDigraph));
	this->currentDigraph = new Digraph(*(donorChain.currentDigraph));

	// Copied manually by iterating through elements.
	// Rewind the chain manually.
	ChainElement* originalList = donorChain.currentElement;
	while(originalList->previous != 0) {
		originalList = originalList->previous;
	}
	// Copy the links until completed
	if (originalList->isCrimp()) {
		this->currentElement = new ChainCrimp(*(ChainCrimp*)originalList);
	} else {
		std::cerr << "The chain you are attempting to copy is invalid -- it does not start with a crimp!" <<std::endl;
	}
	while(originalList->next != 0) {
		originalList = originalList->next;
		// Select the correct copy constructor based on what kind of element it is:
		if (originalList->isCrimp()) {
			this->currentElement->next = new ChainCrimp(*(ChainCrimp*)originalList);
		} else {
			this->currentElement->next = new ChainLink(*(ChainLink*)originalList);
		}
		this->currentElement->next->previous = this->currentElement;
		this->currentElement = this->currentElement->next;
	}
	// Be kind, rewind
	this->currentElementIndex = 0;
	this->resetX();
	while(this->currentElement->previous != 0) {
		this->currentElement = this->currentElement->previous;
	}
}

Chain::~Chain() {
	delete(this->startingDigraph);
	this->clearChain();
	delete(this->currentElement);
	delete(this->linkEvaluationFunction);
	delete(this->crimpEvaluationFunction);
	delete(this->currentDigraph);
	delete(this->betaScore);
	delete(this->betaInformation);
	#if GAMMA_LENGTH > 0
	delete(this->gammaScore);
	delete(this->gammaInformation);
	#endif
}

void Chain::insertCrimpAfter(Digraph* observedDigraph) {
	ChainCrimp *newCrimp = new ChainCrimp(this->actors);
	newCrimp->observedDigraph = observedDigraph;
	if (this->currentElement) {
		newCrimp->previous = this->currentElement;
		ChainElement* temporary = this->currentElement->next;
		this->currentElement->next = newCrimp;
		newCrimp->next = temporary;
		this->advanceChain();
	} else {
		this->currentElement = newCrimp;
	}
	this->chainSize++;
	this->crimps++;
}

void Chain::insertCrimpAfter(const char* fileName) {
	Digraph *newDigraph = new Digraph(this->actors);
	newDigraph->readNet(fileName);
	this->insertCrimpAfter(newDigraph);
}

void Chain::insertCrimpAfter() {
	Digraph *newDigraph = new Digraph(*(this->currentDigraph));
	this->insertCrimpAfter(newDigraph);
}

void Chain::evaluateChain() {
	this->evaluateChain(0,this->chainSize-1);
}

void Chain::evaluateChain(int start, int stop) {
	if (start > stop || start < 0 || stop >= this->chainSize) {
		std::cerr << "Error in evaluate chain from " << start << " -- " << stop << std::endl;
	}
	this->resetLikelihoodStatistics();
	bool flag = false;
	if (this->currentElement != 0) {
		this->rewind();
		flag = true;
	}
	while (flag) {
		if ((this->currentElementIndex>=start) | (this->currentElementIndex <=stop)) {
			if (this->currentElement->isLink()) {
				this->currentDigraph->flip( ((ChainLink*)this->currentElement)->i, ((ChainLink*)this->currentElement)->j);
				this->currentElement->evaluate(this->currentDigraph, this->linkEvaluationFunction);
				this->currentDigraph->flip( ((ChainLink*)this->currentElement)->i, ((ChainLink*)this->currentElement)->j);
				this->incrementLogLikelihood(this->currentElement->getLogProbability());
				for (int i=0; i<BETA_LENGTH; i++) {
					this->incrementBetaScore(i, this->currentElement->getScore(i));
					for (int j=0; j<BETA_LENGTH; j++) {
						this->incrementBetaInformation(i, j, this->currentElement->getInformation(i,j) );
					}
				}
			} else {
				this->currentElement->evaluate(this->currentDigraph, this->crimpEvaluationFunction);
				this->incrementLogLikelihood(this->currentElement->getLogProbability());
				for (int i=0; i<GAMMA_LENGTH; i++) {
					this->incrementGammaScore(i, this->currentElement->getScore(i));
					for (int j=0; j<GAMMA_LENGTH; j++) {
						this->incrementGammaInformation(i, j, this->currentElement->getInformation(i,j) );
					}
				}
			}
		} else {
			std::cerr << "Error in evaluateChain: the currentElementIndex is out of bounds!" << std::endl;
		}
		if (this->currentElement->next) {
			this->advanceChain();
		} else {
			flag = false;
		}
	}
}

void Chain::resetLikelihoodStatistics() {
	this->setLogLikelihood(0);
	for (int i=0; i<BETA_LENGTH; i++) {
		this->setBetaScore(i, 0);
		for (int j=0; j<BETA_LENGTH; j++) {
			this->setBetaInformation(i,j,0);
		}
	}
	for (int i=0; i<GAMMA_LENGTH; i++) {
		this->setGammaScore(i, 0);
		for (int j=0; j<GAMMA_LENGTH; j++) {
			this->setGammaInformation(i,j,0);
		}
	}
}

void Chain::resetChain() {
	this->clearChain();
	this->resetX();
	this->resetLikelihoodStatistics();
}

void Chain::resetStart() {
	for (int i=0; i<this->actors; i++) {
		for (int j=0; j<this->actors; j++) {
			this->startingDigraph->set(i,j,false);
		}
	}
}

int Chain::getChainSize() {
	return this->chainSize;
}

void Chain::clearChain() {
	this->fastForward();
	while (this->currentElement != 0) {
		ChainElement* temporary = this->currentElement;
		this->retreatChain();
		delete(temporary);
	}
	this->diagonalElements=0;
	this->currentElement = 0;
	this->currentElementIndex=0;
	this->chainSize=0;
	this->crimps=0;
}

void Chain::resetX() {
	this->rewind();
	for (int i=0; i<this->actors; i++) {
		for (int j=0; j<this->actors; j++) {
			this->currentDigraph->set(i,j, this->startingDigraph->get(i,j));
		}
	}
}

void Chain::simulateStartingDigraph(int steps, MersenneTwister *r) {
	ChainLink s(this->actors);
	while (steps-- >0) {
		s.simulate(this->startingDigraph,this->linkEvaluationFunction,r);
		this->startingDigraph->flip(s.i,s.j);
	}
}

int Chain::getActors() {
	return this->actors;
}

int Chain::getDiagonalElements() {
	return this->diagonalElements;
}

void Chain::rewind() {
	if (this->currentElement != 0) {
		while (this->currentElement->previous != 0) {
			this->retreatChain();
		}
	}
}

void Chain::fastForward() {
	if (this->currentElement != 0) {
		while (this->currentElement->next != 0) {
			this->advanceChain();
		}
	}
}

void Chain::setPosition(int position) {
	if (this->currentElement != 0) {
		if (position > this->chainSize || position < 0) {
			std::cerr << "You requested position " << position
					<< " of a chain of size " << this->chainSize << "!" << std::endl;
		} else {
			while ((position > currentElementIndex) & (this->currentElement->next != 0)) {
				this->advanceChain();
			}
			while ((position < currentElementIndex) & (this->currentElement->previous != 0)) {
				this->retreatChain();
			}
		}
		if (position != this->currentElementIndex) {
			std::cerr << "Something went wrong with setting position "
					<< position << ", you are at " << this->currentElementIndex << std::endl;
		}
	}
}

void Chain::insertLinkAfter(int i, int j) {
	if (this->currentElement) {
		ChainLink* appendedLink = new ChainLink(this->actors);
		appendedLink->i = i;
		appendedLink->j = j;
		appendedLink->previous = this->currentElement;
		if (this->currentElement->next) {
			appendedLink->next = this->currentElement->next;
			currentElement->next->previous = appendedLink;
			this->currentElement->next = appendedLink;
		} else {
			this->currentElement->next = appendedLink;
		}
		this->chainSize++;
		if ((i == j) & (i >= 0)) {
			this->diagonalElements++;
		}
		this->currentElement->next->evaluate(this->currentDigraph, this->linkEvaluationFunction);
		this->advanceChain();
		for (int i=0; i<BETA_LENGTH; i++) {
			this->incrementBetaScore(i, this->currentElement->getScore(i));
			for (int j=0; j<BETA_LENGTH; j++) {
				this->incrementBetaInformation(i,j,this->currentElement->getInformation(i,j));
			}
		}
		this->incrementLogLikelihood(this->currentElement->getLogProbability());
	} else {
		std::cerr << "You cannot insert a link until you have added a crimp to the chain." << std::endl;
	}
}

void Chain::deleteCurrentElement() {
	if(this->currentElement == 0) {
		std::cerr << "Tried to delete a null element!" << std::endl;
	} else if ( this->currentElement->isCrimp() & this->currentElementIndex == 0) {
		std::cerr << "Error: You are deleting the first crimp!" << std::endl;
	} else if (this->currentElementIndex == 0) {
		std::cerr << "Error: The first element is NOT a crimp, AND you tried to delete it!" << std::endl;
	} else {
		if (this->currentElement->isLink()) {
			if (((ChainLink*)this->currentElement)->i == ((ChainLink*)this->currentElement)->j) {
				this->diagonalElements--;
			}
		} else if (this->currentElement->isCrimp()) {
			this->crimps--;
		}
		if (this->currentElement->previous) {
			this->currentElement->previous->next = this->currentElement->next;
		}
		if (this->currentElement->next) {
			this->currentElement->next->previous = this->currentElement->previous;
		}
		ChainElement* temp = 0;
		if (this->currentElement->previous) {
			 temp = this->currentElement->previous;
			this->currentElementIndex--;
		} else if (this->currentElement->next) {
			temp = this->currentElement->next;
		}
		delete(this->currentElement);
		this->currentElement = temp;
		this->chainSize--;
	}
}

void Chain::incrementBetaScore(int a, double value) {
	this->setBetaScore(a, this->getBetaScore(a) + value );
}

void Chain::incrementBetaInformation(int a, int b, double value) {
	if ((a<BETA_LENGTH)&(b<BETA_LENGTH)) {
		this->setBetaInformation(a,b, this->getBetaInformation(a, b) + value );
	} else {
		std::cerr << "Chain::setBetaInformation Error: index out of bounds " << a << "," << b << std::endl;
	}
}

void Chain::setBetaScore(int a, double value) {
	if ((a<BETA_LENGTH)) {
		this->betaScore[a] = value;
	} else {
		std::cerr << "Chain::setBetaScore Error: index out of bounds " << a << std::endl;
	}
}

void Chain::setBetaInformation(int a, int b, double value) {
	if ((a<BETA_LENGTH)&(b<BETA_LENGTH)) {
		this->betaInformation[a * BETA_LENGTH + b] = value;
	} else {
		std::cerr << "Chain::setBetaInformation Error: index out of bounds " << a << "," << b << std::endl;
	}
}

void Chain::printLinkedList() {
	std::cout << "Printing linked list references:" << std::endl;
	this->rewind();
	if (this->currentElement) {
		while (this->currentElement->next != 0) {
			std::cout << this->currentElementIndex << ") " <<
					this->currentElement->previous << " -> " <<
					this->currentElement << " -> " <<
					this->currentElement->next << std::endl;
			this->advanceChain();
		}
		std::cout << this->currentElementIndex << ") " <<
				this->currentElement->previous << " -> " <<
				this->currentElement << " -> " <<
				this->currentElement->next << std::endl;
	} else {
		std::cout << "  *Chain is empty." <<std::endl;
	}
}

void Chain::print() {
	int thisPosition = currentElementIndex;
	std::cout << "Starting digraph:" << std::endl;
	this->startingDigraph->print();
	this->rewind();
	std::cout << "Chain: " << this->chainSize << " elements, " << this->diagonalElements << " diagonal,"
			<< this->crimps << " observations." << std::endl;
	std::cout << " (IND)   i   j P:i->j  | {Scores|Information|Statistics} " << std::endl;
	if (this->currentElement != 0) {
		while (this->currentElement->next != 0) {
			this->currentElement->print(this->currentElementIndex);
			this->advanceChain();
		}
		this->currentElement->print(this->currentElementIndex);
	}
	std::cout << "Ending digraph:" << std::endl;
	this->currentDigraph->print();
	this->setPosition(thisPosition);
	printf("LogLike[(chain)]=%5.0f; Score=", this->logLikelihood);
	for (int i=0; i<BETA_LENGTH; i++) {
		printf(" %+0.3f",this->getBetaScore(i) );
	}
	printf("; Info=");
	for (int i=0; i<BETA_LENGTH; i++) {
		printf(" %+0.3f", this->getBetaInformation(i,i));
	}
	printf("\n");

}

void Chain::advanceChain() {
	if (this->currentElement) {
		this->currentElementIndex++;
		this->currentElement = this->currentElement->next;
		if (this->currentElement->isLink()) {
			this->currentDigraph->flip(((ChainLink*)this->currentElement)->i, ((ChainLink*)this->currentElement)->j);
		} else {
			// The crimp does not modify the true digraph
		}
	}
}

void Chain::retreatChain() {
	if (this->currentElement) {
		if (this->currentElement->isLink()) {
			this->currentDigraph->flip(((ChainLink*)this->currentElement)->i, ((ChainLink*)this->currentElement)->j);
		} else {
			// The crimp does not modify the true digraph
		}
		this->currentElementIndex--;
		this->currentElement = this->currentElement->previous;
	}
}

void Chain::insertSimulatedLinkAfter(MersenneTwister* r) {
	ChainLink temporary(this->actors);
	temporary.simulate(this->currentDigraph, this->linkEvaluationFunction, r);
	this->insertLinkAfter(temporary.i, temporary.j);
}

void Chain::insertSimulatedCrimpAfter(MersenneTwister* r) {
	Digraph* observedDigraph = new Digraph(this->actors);
	ChainCrimp temporary(this->actors);
	temporary.simulate(this->currentDigraph, observedDigraph, this->crimpEvaluationFunction, r);
	this->insertCrimpAfter(observedDigraph);
}

void Chain::incrementGammaScore(int a, double value) {
	this->setGammaScore(a, this->getGammaScore(a) + value);
}

void Chain::incrementGammaInformation(int a, int b, double value) {
	this->setGammaInformation(a,b,this->getGammaInformation(a,b) + value);
}

void Chain::setGammaScore(int a, double value) {
	if ((a<GAMMA_LENGTH)) {
		this->gammaScore[a] = value;
	} else {
		std::cerr << "Chain::getGammaScore Error: index out of bounds " << a << std::endl;
	}
}

double Chain::getScore(int a) {
	if (a < BETA_LENGTH) {
		return this->getBetaScore(a);
	} else {
		return this->getGammaScore(a-BETA_LENGTH);
	}
}

double Chain::getInformation(int a, int b) {
	if ((a < BETA_LENGTH) & (b < BETA_LENGTH)) {
		return this->getBetaInformation(a, b);
	} else if ((a >= BETA_LENGTH) & (b >= BETA_LENGTH) & (a < THETA_LENGTH) & (b < THETA_LENGTH)){
		return this->getGammaInformation(a-BETA_LENGTH, b-BETA_LENGTH);
	} else {
		return 0;
	}
}

void Chain::setGammaInformation(int a, int b, double value) {
	if ((a<GAMMA_LENGTH)&(b<GAMMA_LENGTH)) {
		this->gammaInformation[a*GAMMA_LENGTH+b] = value;
	} else {
		std::cerr << "Chain::setGammaInformation Error: index out of bounds " << a << "," << b << std::endl;
	}
}

double Chain::getGammaScore(int k) {
	if ((k<GAMMA_LENGTH)) {
		return this->gammaScore[k];
	} else {
		std::cerr << "Chain::getGammaScore Error: index out of bounds " << k << std::endl;
		return 0.0;
	}
}

double Chain::getGammaInformation(int a, int b) {
	if ((a<GAMMA_LENGTH)&(b<GAMMA_LENGTH)) {
		return this->gammaInformation[a*GAMMA_LENGTH + b];
	} else {
		std::cerr << "Chain::getGammaInformation Error: index out of bounds " << a << "," << b << std::endl;
		return 0.0;
	}
}

double Chain::getBetaScore(int k) {
	if ((k<BETA_LENGTH)) {
		return this->betaScore[k];
	} else {
		std::cerr << "Chain::getBetaScore Error: index out of bounds " << k << std::endl;
		return 0.0;
	}
}

double Chain::getBetaInformation(int a, int b) {
	if ((a<BETA_LENGTH)&(b<BETA_LENGTH)) {
		return this->betaInformation[a * BETA_LENGTH + b];
	} else {
		std::cerr << "Chain::getBetaInformation Error: index out of bounds " << a << "," << b << std::endl;
		return 0.0;
	}
}

double Chain::getLogLikelihood() {
	return this->logLikelihood;
}

void Chain::setLogLikelihood(double value) {
	this->logLikelihood = value;
}

void Chain::incrementLogLikelihood(double value) {
	this->logLikelihood += value;
}

bool Chain::drawProposal(MersenneTwister* r) {
	double selector = r->next();
	bool result = false;

	if ((selector -= PERMUTE_PROPOSAL) < 0) {
		// Permute
		result = this->permute(r);
	} else if ((selector -= INSERTD_PROPOSAL) < 0) {
		// Insert diagonal
		result = this->insertDiagonal(r);
	} else if ((selector -= DELETED_PROPOSAL) < 0) {
		// Delete diagonal
		result = this->deleteDiagonal(r);
	} else if ((selector -= INSERTCCP_PROPOSAL) < 0) {
		// Insert Consecutively Cancelling Pairs
		result = this->insertConsecutivelyCancelingPair(r);
	} else if ((selector -= DELETECCP_PROPOSAL) < 0) {
		// Delete Consecutively Cancelling Pairs
		result = this->deleteConsecutivelyCancelingPair(r);
	} else if ((selector -= INSERTEC_PROPOSAL) < 0) {
		// Delete Consecutively Cancelling Pairs
		result = this->insertEndCrimpElement(r);
	} else if ((selector -= DELETEEC_PROPOSAL) < 0) {
		// Delete Consecutively Cancelling Pairs
		result = this->deleteEndCrimpElement(r);
	} else {
		std::cerr << "Your selection probabilities for proposal type do not add to one." << std::endl;
	}
	return result;
}

bool Chain::insertDiagonal(MersenneTwister* r) {
	// Decide on a position to insert the diagonal into
	int insertionIndex = (this->chainSize-2) * r->next() + 1;
	this->setPosition(insertionIndex);
	int ego = this->actors * r->next();
	ChainLink proposedLink(this->actors);
	proposedLink.i = ego;
	proposedLink.j = ego;
	double logProb = proposedLink.evaluateProposal(this->currentDigraph,this->linkEvaluationFunction);

	double acceptanceProbability = this->getAlpha()  * exp(logProb)
			* (DELETED_PROPOSAL / INSERTD_PROPOSAL) / (double)(this->diagonalElements + 1) ;
	if ( r->next() < acceptanceProbability ) {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Accepting proposal to insert element (" << ego << "," << ego << ") at position "
				<< this->currentElementIndex + 1 << " p=" << acceptanceProbability << std::endl;
#endif
		this->insertLinkAfter(ego,ego);
		return true;
	} else {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Rejecting proposal to insert element (" << ego << "," << ego << ") at position "
				<< this->currentElementIndex + 1 << " p=" << acceptanceProbability << std::endl;
#endif
		return false;
	}
}

int Chain::getLinkCount() {
	return this->chainSize - this->crimps;
}

bool Chain::deleteDiagonal(MersenneTwister* r) {
	if (this->diagonalElements > 0) {
		int deletionIndex = this->diagonalElements * r->next();
		this->setPosition(0);
		while (deletionIndex>=0) {
			this->advanceChain();
			if (this->currentElement->isLink()) {
				if (((ChainLink*)this->currentElement)->i == ((ChainLink*)this->currentElement)->j) {
					deletionIndex--;
				}
			}
		}
		if (this->currentElement->isCrimp()) {
			std::cerr << "Something is wrong with the delete diagonal proposal -- ended at a crimp!" << std::endl;
		} else if (((ChainLink*)this->currentElement)->i != ((ChainLink*)this->currentElement)->j) {
			std::cerr << "Something is wrong with the delete diagonal proposal -- this link is not diagonal!" << std::endl;
		}
		double acceptanceProbability = (1/this->getAlpha()) * exp(-1 * this->currentElement->getLogProbability() ) *
				(double)(this->diagonalElements) ;
		if ( r->next() < acceptanceProbability ) {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Accepting proposal to delete element (" << ((ChainLink*)this->currentElement)->i
				<< ","<< ((ChainLink*)this->currentElement)->i << ") at "
				<< this->currentElementIndex << " p=" << acceptanceProbability << std::endl;
#endif
			this->deleteCurrentElement();
			return true;
		} else {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Rejecting proposal to delete element (" << ((ChainLink*)this->currentElement)->i << ","
				<< ((ChainLink*)this->currentElement)->i << ") at "
				<< this->currentElementIndex << " p=" << acceptanceProbability << std::endl;
#endif
			return false;
		}
	} else {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Rejecting proposal to delete element, as there are no diagonal elements to delete." << std::endl;
#endif
		return false;
	}
}

bool Chain::permute(MersenneTwister* r) {
	// Decide on the permutation size, ensures the interval [2,PERMUTE_MAX]
	int permutationSize = (PERMUTE_MAX - 1) * r->next() + 2;
	// Check for a really short chain, and adjust permutation size.
	if (permutationSize > this->chainSize - 2) {
		permutationSize = this->chainSize - 2;
		std::cerr << "Chain was too short for the proposed size. Please scale PERMUTE_MAX back--this will affect MCMC convergence!" << std::endl;
	}
	// Determine a starting value on [1, (chainSize - 1) - (permutationSize-1)]
	int start = (this->chainSize - permutationSize - 2) * r->next() + 1;
	int end = start + permutationSize-1;
	if ((start==0) | (end==this->chainSize-1) | ((end-start) > permutationSize) | (start==end)) {
		std::cerr << "Something is wrong with the coding of Chain::permute--start and end were not properly generated: (" << start << "," << ")" << std::endl;
	}
	// Determine the log odds for the portion of the chain under consideration:
	this->setPosition(start);
	double originalLogOdds = 0;
	int crimpsIncluded = 0;
	for (int i=0; i+start <= end; i++) {
		this->permuteProposal[i] = this->currentElement;
		originalLogOdds+= this->currentElement->getLogProbability();
		if (this->currentElement->isCrimp()) {
			crimpsIncluded++;
		}
		if (crimpsIncluded > 1) {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Rejecting proposal to permute " << start << " to " << end << ", p=NA, included 2 crimps!" << std::endl;
#endif
			return false;
		}
#if GAMMA_LENGTH == 0
		if (crimpsIncluded > 0) {
	#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Rejecting proposal to permute " << start << " to " << end << ", p=NA, included 1 crimps and no IA model included!" << std::endl;
	#endif
			return false;
		}
#endif
		this->advanceChain();
	}
	// Decide on a permutation of the indices
	int p[permutationSize];
	this->permutedIntegers(permutationSize, r, p);
	double proposedLogOdds = 0;
	// Need to set the permutation to before the start point
	this->setPosition(start-1);
	// std::cout << "Proposed permutation: " << std::endl;
	for (int i=0; i+start <= end; i++) {
		if (this->permuteProposal[p[i]]->isLink()) {
			((ChainLink*)this->permuteProposal[p[i]])->evaluate(this->currentDigraph, this->linkEvaluationFunction);
			this->currentDigraph->flip(((ChainLink*)this->permuteProposal[p[i]])->i, ((ChainLink*)this->permuteProposal[p[i]])->j);
		} else if (this->permuteProposal[p[i]]->isCrimp()) {
			((ChainCrimp*)this->permuteProposal[p[i]])->evaluate(this->currentDigraph, this->crimpEvaluationFunction );
			// Crimp does not modify the digraph
		}
		proposedLogOdds += this->permuteProposal[p[i]]->getLogProbability();
	}
	// Reset X to the start position
	for (int i=0; i+start <= end; i++) {
		if (this->permuteProposal[p[i]]->isLink()) {
			this->currentDigraph->flip(((ChainLink*)this->permuteProposal[p[i]])->i, ((ChainLink*)this->permuteProposal[p[i]])->j);
		} else {
			// Crimp does not modify the digraph
		}
	}
	double acceptanceProbability = exp(proposedLogOdds - originalLogOdds);
	if (r->next() < acceptanceProbability ) {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Accepting proposal to permute " << start << " to " << end << ", p=" << acceptanceProbability << std::endl;
#endif
		// Swap the proposals for the current elements, i.e. shuffle the elements in the chain.
		this->setPosition(start-1);
		ChainElement* pre = this->currentElement;
		this->setPosition(end+1);
		ChainElement* post = this->currentElement;
		pre->next = this->permuteProposal[p[0]];
		this->permuteProposal[p[0]]->previous = pre;
		this->permuteProposal[p[0]]->next = this->permuteProposal[p[1]];
		for (int i=1; i < end - start; i++) {
			// Make the permutations internally consistent
			this->permuteProposal[p[i]]->previous = this->permuteProposal[p[i-1]];
			this->permuteProposal[p[i]]->next = this->permuteProposal[p[i+1]];
		}

		this->permuteProposal[p[end-start]]->previous = this->permuteProposal[p[end-start-1]];
		this->permuteProposal[p[end-start]]->next = post;
		post->previous = this->permuteProposal[p[end-start]];
		//this->evaluateChain(start,end);
		return true;
	} else {
#ifdef PRINT_PROPOSAL_RESULTS
		std::cout << "**Rejecting proposal to permute " << start << " to " << end << ", p=" << acceptanceProbability << std::endl;
#endif
		return false;
	}
}

bool Chain::insertConsecutivelyCancelingPair(MersenneTwister* r) {
	//TODO: Implement insertConsecutivelyCancelingPair.
	return false;
}

bool Chain::deleteConsecutivelyCancelingPair(MersenneTwister* r) {
	//TODO: Implement deleteConsecutivelyCancelingPair.
	return false;
}

bool Chain::insertEndCrimpElement(MersenneTwister* r) {
	//TODO: Implement insertEndCrimpElement.
	return false;
}

bool Chain::deleteEndCrimpElement(MersenneTwister* r) {
	//TODO: Implement deleteEndCrimpElement.
	return false;
}

void Chain::permutedIntegers (int size, MersenneTwister* r, int *integers) {
	double x[size];
	// Generate a random vector of doubles
	for (int i=0; i<size; i++) {
		x[i] = r->next();
		integers[i] = 0;
	}
	// Sort the vector, and adjust the integer vectors
	for (int i=0; i<size; i++) {
		for (int j = 0; j<size; j++) {
			if (x[i] > x[j]) {
				integers[i]++;
	// Only happens once in a blue moon, due to rounding
			} else if ((x[i] == x[j]) & (i != j)) {
				x[i]+=.000000001;
				integers[i]++;
			}
		}
	}
}

void Chain::incrementAlpha(double value) {
	this->setAlpha(this->alpha + value);
}

void Chain::setAlpha(double value) {
	if (value > 0) {
		this->alpha = value;
	} else {
		std::cout << "Chain::setAlpha: Cannot have basic rate parameter less than zero:" << value << std::endl;
	}
}

double Chain::getAlpha() {
	return this->alpha;
}

int Chain::getCrimpCount() {
	return this->crimps;
}

double Chain::getAlphaScore() {
	return ((double)this->chainSize / this->alpha - (double)this->actors);
}

double Chain::getAlphaInformation() {
	return ((double)this->chainSize / (this->alpha * this->alpha));
}

}
